
#######################################
#### ROBUSTNESS CHECKS FOR TABLE 4 ####
#######################################

# library(pacman)
# p_load(dplyr, readr, fixest, modelsummary, plm, kableExtra, stevemisc, lmtest) 
library(dplyr)
library(readr)
library(fixest)
library(modelsummary)
library(plm)
library(kableExtra)
library(stevemisc)
library(lmtest)

options(modelsummary_factory_default = 'kableExtra')

alldata1 <- read_rds("replication_data1.rds")
alldata2 <- read_rds("replication_data2.rds")
alldata3 <- read_rds("replication_data3.rds")

#### F9: ALL CONTROL VARIABLES ####

# In the paper, only the main independent variable is shown. 
cm <- c("weight" = "ln(TC connections)",
        "patents_cty_1" = "Patents (exporter) as share of GDP",
        "patents_cty_2" = "Patents (importer) as share of GDP",
        "gdpcap_d" = "GDP per capita (exporter)", "gdpcap_o" = "GDP per capita (importer)",
        "log(pop_d)" = "Ln(population) (exporter)", "log(pop_o)" = "Ln(population) (importer)",
        "distcap" = "Distance between capitals", "contig" = "Contiguity",
        "comlang_off" = "Common language", "fta_wto" = "Regional trade agreement", "wto_dyad" = "WTO dyad",
        "dem_dyad" = "Democratic dyad", "pta" = "Preferential trade agreement", "comcur" = "Common currency",
        "alliance" = "Alliance",
        "rivalry" = "Strategic rivals")

mod_all0 <- feols(logtrade ~ weight | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1)

mod_all2 <- feols(logtrade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2)

mod_all4 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex", # This line produces latex output of the table
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects by dyad, country and year, clustered standard errors by dyad and year.",
                                     "Zero imputation on dyads with missing on TC connections.",
                                     "Gravity controls: GDP, population, distance between capitals, common language, regional trade agreement, WTO dyad.",
                                     "Gravity+R&D controls: Adds to Gravity patents per country.",
                                     "Gravity+ controls: Adds to Gravity+ democratic dyad, preferential trade agreement, common currency.",
                                     "Gravity++ controls: Adds to Gravity++ strategic rivalry, alliance."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 8, full_width = F) 

#### F10: IMF DEP VAR ####

cm <- c("weight" = "ln(TC connections)", 
        "patents_cty_1" = "Patents (exporter) as share of GDP",
        "patents_cty_2" = "Patents (importer) as share of GDP",
        "gdpcap_d" = "GDP per capita (exporter)", "gdpcap_o" = "GDP per capita (importer)",
        "log(pop_d)" = "Ln(population) (exporter)", "log(pop_o)" = "Ln(population) (importer)",
        "distcap" = "Distance between capitals", "contig" = "Contiguity",
        "comlang_off" = "Common language", "fta_wto" = "Regional trade agreement", "wto_dyad" = "WTO dyad",
        "dem_dyad" = "Democratic dyad", "pta" = "Preferential trade agreement", "comcur" = "Common currency",
        "alliance" = "Alliance",
        "rivalry" = "Strategic rivals")

mod_all0 <- feols(log_imf_trade ~ weight | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all1 <- feols(log_imf_trade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1)

mod_all2 <- feols(log_imf_trade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all3 <- feols(log_imf_trade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2)

mod_all4 <- feols(log_imf_trade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           statistic = c("{std.error}"),
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects: Dyad, country and year.", 
                                     "Clustered standard errors by dyad and year.",
                                     "Zero imputation on dyads with missing on TC connections."),
                           add_rows = rows,
                           coef_omit = c("Intercept|factor(year)*|factor(cty_1)*|factor(cty_2)*|factor(dyad_id)*"))

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (IMF)." = 5)) %>%
  kable_styling(font_size = 10, full_width = FALSE) 

#### F11: SHARE OF TRADE ####

mod_all0 <- feols(share_trade ~ weight | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all1 <- feols(share_trade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1)

mod_all2 <- feols(share_trade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all3 <- feols(share_trade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2)

mod_all4 <- feols(share_trade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 5,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           statistic = c("{std.error}"),
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects: Dyad and year.", 
                                     "Clustered standard errors by dyad and year.",
                                     "Zero imputation on dyads with missing on TC connections."),
                           add_rows = rows,
                           coef_omit = c("Intercept|factor(year)*|factor(cty_1)*|factor(cty_2)*|factor(dyad_id)*"))

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: Share of trade" = 5)) %>%
  kable_styling(font_size = 10, full_width = FALSE)

#### F12: DICHOTOMOUS DEPENDENT VARIABLE ####

alldata1_dich <- alldata1 %>%
  mutate(weight = ifelse(weight >= 1, 1, weight))

alldata2_dich <- alldata2 %>%
  mutate(weight = ifelse(weight >= 1, 1, weight))

alldata3_dich <- alldata3 %>%
  mutate(weight = ifelse(weight >= 1, 1, weight))

mod_all0 <- feols(logtrade ~ weight | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1_dich)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1_dich)

mod_all2 <- feols(logtrade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1_dich)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2_dich)

mod_all4 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3_dich)

cm <- c("weight" = "Presence of TC connection",
        "patents_cty_1" = "Patents (exporter) as share of GDP",
        "patents_cty_2" = "Patents (importer) as share of GDP",
        "gdpcap_d" = "GDP per capita (exporter)", "gdpcap_o" = "GDP per capita (importer)",
        "log(pop_d)" = "Ln(population) (exporter)", "log(pop_o)" = "Ln(population) (importer)",
        "distcap" = "Distance between capitals", "contig" = "Contiguity",
        "comlang_off" = "Common language", "fta_wto" = "Regional trade agreement", "wto_dyad" = "WTO dyad",
        "dem_dyad" = "Democratic dyad", "pta" = "Preferential trade agreement", "comcur" = "Common currency",
        "alliance" = "Alliance",
        "rivalry" = "Strategic rivals")

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects: Dyad and year.", 
                                     "Clustered standard errors by dyad and year."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: Presence of TC membership (binary) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 10, full_width = FALSE) 

#### F13: WITH COUNTRY FIXED EFFECTS ####

cm <- c("weight" = "ln(TC connections)",
        "patents_cty_1" = "Patents (exporter) as share of GDP",
        "patents_cty_2" = "Patents (importer) as share of GDP",
        "gdpcap_d" = "GDP per capita (exporter)", "gdpcap_o" = "GDP per capita (importer)",
        "log(pop_d)" = "Ln(population) (exporter)", "log(pop_o)" = "Ln(population) (importer)",
        "distcap" = "Distance between capitals", "contig" = "Contiguity",
        "comlang_off" = "Common language", "fta_wto" = "Regional trade agreement", "wto_dyad" = "WTO dyad",
        "dem_dyad" = "Democratic dyad", "pta" = "Preferential trade agreement", "comcur" = "Common currency",
        "alliance" = "Alliance",
        "rivalry" = "Strategic rivals")

mod_all0 <- feols(logtrade ~ weight | 
                    cty_1+cty_2+year, 
                  cluster = "cty_1+cty_2+year", data = alldata1)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    cty_1+cty_2+year,
                  cluster = "cty_1+cty_2+year", data = alldata1)

mod_all2 <- feols(logtrade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    cty_1+cty_2+year, 
                  cluster = "cty_1+cty_2+year", data = alldata1)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    cty_1+cty_2+year, 
                  cluster = "cty_1+cty_2+year", data = alldata2)

mod_all4 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    cty_1+cty_2+year, 
                  cluster = "cty_1+cty_2+year", data = alldata3)


models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects: Country 1, country 2 and year.", 
                                     "Clustered standard errors by country 1, country 2 and year.",
                                     "Zero imputation on dyads with missing on TC connections."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 8, full_width = FALSE) 

#### F14: WITH REGION FIXED EFFECTS ####

mod_all0 <- feols(logtrade ~ weight | 
                    cty_1_region+cty_2_region+year, 
                  cluster = "cty_1_region+cty_2_region+year", data = alldata1)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    cty_1_region+cty_2_region+year,
                  cluster = "cty_1_region+cty_2_region+year", data = alldata1)

mod_all2 <- feols(logtrade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    cty_1_region+cty_2_region+year, 
                  cluster = "cty_1_region+cty_2_region+year", data = alldata1)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    cty_1_region+cty_2_region+year, 
                  cluster = "cty_1_region+cty_2_region+year", data = alldata2)

mod_all4 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    cty_1_region+cty_2_region+year, 
                  cluster = "cty_1_region+cty_2_region+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects: Region for country 1, region for country 2 and year.", 
                                     "Clustered standard errors by region for country 1, region for country 2 and year.",
                                     "Zero imputation on dyads with missing on TC connections."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 10, full_width = FALSE) 

#### F15: WITH ROLLING TIME TREND ####

alldata1_rollmean <- alldata1 %>%
  mutate(period = cut(alldata1$year, seq(2000, 2025, by = 5), right = FALSE)) %>% 
  dplyr::group_by(dyad_id, period) %>% 
  dplyr::summarise(logtrade = mean(logtrade), 
                   weight = mean(weight),
                   gdpcap_d =mean(gdpcap_d),
                   gdpcap_o = mean(gdpcap_o),
                   pop_d = mean(pop_d),
                   pop_o = mean(pop_o),
                   distcap = mean(distcap),
                   comcol = mean(comcol),
                   comlang_off = mean(comlang_off),
                   fta_wto = mean(fta_wto),
                   wto_dyad = mean(wto_dyad),
                   patents_cty_1 = mean(patents_cty_1),
                   patents_cty_2 = mean(patents_cty_2)) %>%
  ungroup()

alldata2_rollmean <- alldata2 %>%
  mutate(period = cut(alldata2$year, seq(2000, 2025, by = 5), right = FALSE)) %>% 
  dplyr::group_by(dyad_id, period) %>% 
  dplyr::summarise(logtrade = mean(logtrade), 
                   weight = mean(weight),
                   gdpcap_d =mean(gdpcap_d),
                   gdpcap_o = mean(gdpcap_o),
                   pop_d = mean(pop_d),
                   pop_o = mean(pop_o),
                   distcap = mean(distcap),
                   comcol = mean(comcol),
                   comlang_off = mean(comlang_off),
                   fta_wto = mean(fta_wto),
                   wto_dyad = mean(wto_dyad),
                   patents_cty_1 = mean(patents_cty_1),
                   patents_cty_2 = mean(patents_cty_2),
                   dem_dyad = mean(dem_dyad),
                   pta = mean(pta),
                   comcur = mean(comcur)) %>%
  ungroup()

alldata3_rollmean <- alldata3 %>%
  mutate(period = cut(alldata3$year, seq(2000, 2025, by = 5), right = FALSE)) %>% 
  dplyr::group_by(dyad_id, period) %>% 
  dplyr::summarise(logtrade = mean(logtrade), 
                   weight = mean(weight),
                   gdpcap_d =mean(gdpcap_d),
                   gdpcap_o = mean(gdpcap_o),
                   pop_d = mean(pop_d),
                   pop_o = mean(pop_o),
                   distcap = mean(distcap),
                   comcol = mean(comcol),
                   comlang_off = mean(comlang_off),
                   fta_wto = mean(fta_wto),
                   wto_dyad = mean(wto_dyad),
                   patents_cty_1 = mean(patents_cty_1),
                   patents_cty_2 = mean(patents_cty_2),
                   dem_dyad = mean(dem_dyad),
                   pta = mean(pta),
                   comcur = mean(comcur),
                   rivalry = mean(rivalry),
                   alliance = mean(alliance)) %>%
  ungroup()

mod_all0 <- feols(logtrade ~ weight | 
                    dyad_id+period, 
                  cluster = "dyad_id+period", 
                  data = alldata1_rollmean)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+period, 
                  cluster = "dyad_id+period", 
                  data = alldata1_rollmean)

mod_all2 <- feols(logtrade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+period, 
                  cluster = "dyad_id+period", 
                  data = alldata1_rollmean)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+period, 
                  cluster = "dyad_id+period", 
                  data = alldata2_rollmean)

mod_all4 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry | 
                    dyad_id+period, 
                  cluster = "dyad_id+period",
                  data = alldata3_rollmean)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects: Region for dyad and 5-year period.", 
                                     "Clustered standard errors by dyad and 5-year period.",
                                     "Zero imputation on dyads with missing on TC connections."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 10, full_width = FALSE) 


#### F16: NO ZERO IMPUTATION ON STANDARD DYADS ####

cm <- c("weight_no_impute" = "ln(TC connections)",
        "patents_cty_1" = "Patents (exporter) as share of GDP",
        "patents_cty_2" = "Patents (importer) as share of GDP",
        "gdpcap_d" = "GDP per capita (exporter)", "gdpcap_o" = "GDP per capita (importer)",
        "log(pop_d)" = "Ln(population) (exporter)", "log(pop_o)" = "Ln(population) (importer)",
        "distcap" = "Distance between capitals", "contig" = "Contiguity",
        "comlang_off" = "Common language", "fta_wto" = "Regional trade agreement", "wto_dyad" = "WTO dyad",
        "dem_dyad" = "Democratic dyad", "pta" = "Preferential trade agreement", "comcur" = "Common currency",
        "alliance" = "Alliance",
        "rivalry" = "Strategic rivals")


mod_all0 <- feols(logtrade ~ weight_no_impute | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all1 <- feols(logtrade ~ weight_no_impute + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1)

mod_all2 <- feols(logtrade ~ weight_no_impute + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all3 <- feols(logtrade ~ weight_no_impute + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2)

mod_all4 <- feols(logtrade ~ weight_no_impute + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects: Dyad and year.", 
                                     "Clustered standard errors by dyad and year."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 10, full_width = FALSE) 


#### F17: WITH PATENTS DATA FROM WIPO

cm <- c("weight" = "ln(TC connections)",
        "wipo_patents_cty_1" = "Patents (exporter) as share of GDP",
        "wipo_patents_cty_2" = "Patents (importer) as share of GDP",
        "gdpcap_d" = "GDP per capita (exporter)", "gdpcap_o" = "GDP per capita (importer)",
        "log(pop_d)" = "Ln(population) (exporter)", "log(pop_o)" = "Ln(population) (importer)",
        "distcap" = "Distance between capitals", "contig" = "Contiguity",
        "comlang_off" = "Common language", "fta_wto" = "Regional trade agreement", "wto_dyad" = "WTO dyad",
        "dem_dyad" = "Democratic dyad", "pta" = "Preferential trade agreement", "comcur" = "Common currency",
        "alliance" = "Alliance",
        "rivalry" = "Strategic rivals")

mod_all0 <- feols(logtrade ~ weight | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1)

mod_all2 <- feols(logtrade ~ weight + wipo_patents_cty_1 + wipo_patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2)

mod_all4 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects by dyad, country and year, clustered standard errors by dyad and year.",
                                     "Zero imputation on dyads with missing on TC connections.",
                                     "Gravity controls: GDP, population, distance between capitals, common language, regional trade agreement, WTO dyad.",
                                     "Gravity+R&D controls: Adds to Gravity patents per country.",
                                     "Gravity+ controls: Adds to Gravity+ democratic dyad, preferential trade agreement, common currency.",
                                     "Gravity++ controls: Adds to Gravity++ strategic rivalry, alliance."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 8, full_width = F) 

#### F18: ONLY GRAVITY CONTROLS, SAME TIME SERIES ####

cm <- c("weight" = "ln(TC connections)",
        "patents_cty_1" = "Patents (exporter) as share of GDP",
        "patents_cty_2" = "Patents (importer) as share of GDP",
        "gdpcap_d" = "GDP per capita (exporter)", "gdpcap_o" = "GDP per capita (importer)",
        "log(pop_d)" = "Ln(population) (exporter)", "log(pop_o)" = "Ln(population) (importer)",
        "distcap" = "Distance between capitals", "contig" = "Contiguity",
        "comlang_off" = "Common language", "fta_wto" = "Regional trade agreement", "wto_dyad" = "WTO dyad",
        "dem_dyad" = "Democratic dyad", "pta" = "Preferential trade agreement", "comcur" = "Common currency",
        "alliance" = "Alliance",
        "rivalry" = "Strategic rivals")

mod_all0 <- feols(logtrade ~ weight | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1)

mod_all2 <- feols(logtrade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad  | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2)

mod_all4 <- feols(logtrade ~ weight  + 
                    fta_wto + wto_dyad |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity", "Gravity")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity", "Gravity",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex",
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects by dyad, country and year, clustered standard errors by dyad and year.",
                                     "Zero imputation on dyads with missing on TC connections."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 8, full_width = F) 

#### F19: GMM MODEL ####

### THIS CODE WILL TAKE VERY LONG TO RUN (APPROX. 24 HOURS) 
### uncomment and run if you have the time

# p_data1 <- pdata.frame(alldata1, index = c("dyad_id", "year"))
# 
# p_data2 <- pdata.frame(alldata2 %>% 
#                          select(logtrade, weight, patents_cty_1, patents_cty_2,
#                                 gdpcap_d, gdpcap_o, pop_d, pop_o,
#                                 distcap, comcol, comlang_off, fta_wto, wto_dyad,
#                                 dem_dyad, pta, comcur, cty_1, cty_2, dyad_id, year) %>%
#                          na.omit(),
#                        index = c("dyad_id", "year"))
# 
# p_data3 <- pdata.frame(alldata3 %>% 
#                          select(logtrade, weight, patents_cty_1, patents_cty_2,
#                                 gdpcap_d, gdpcap_o, pop_d, pop_o,
#                                 distcap, comcol, comlang_off, fta_wto, wto_dyad,
#                                 dem_dyad, pta, comcur, rivalry, alliance, cty_1, cty_2, dyad_id, year) %>%
#                          na.omit(),
#                        index = c("dyad_id", "year"))
# 
# gmm_all0 <- pgmm(logtrade ~ lag(logtrade, 1:2) + lag(weight, 0:1) | lag(logtrade, 2:10),
#                  data = p_data1,
#                  effect="twoways", model="twosteps", transformation = "ld")
# 
# gmm_all1 <- pgmm(logtrade ~ lag(logtrade, 1:2) + lag(weight, 0:1) +
#                    fta_wto + wto_dyad | lag(logtrade, 2:10),
#                  data = p_data1,
#                  effect="individual", model="twosteps", transformation = "ld")
# 
# gmm_all2 <- pgmm(logtrade ~ lag(logtrade, 1:2) + lag(weight, 0:1) + 
#                    patents_cty_1 + patents_cty_2 +
#                    fta_wto + wto_dyad | lag(logtrade, 3:10),
#                  data = p_data1,
#                  effect="individual", model="twosteps", transformation = "ld")
# 
# gmm_all3 <- pgmm(logtrade ~ lag(logtrade, 1:2) + lag(weight, 0:1) + 
#                    fta_wto + wto_dyad +
#                    dem_dyad + pta + comcur | lag(logtrade, 3:10),
#                  data = p_data2,
#                  effect="twoways", model="twosteps")
# 
# gmm_all4 <- pgmm(logtrade ~ lag(logtrade, 1:2) + lag(weight, 0:1) + 
#                    fta_wto + wto_dyad +
#                    dem_dyad + pta + comcur +
#                    alliance + rivalry | lag(logtrade, 3:10),
#                  data = p_data3,
#                  effect="twoways", model="twosteps", transformation = "ld")
# 
# # Serial correlation
# mtest(gmm_all0)
# # Focus on AR 2 (second lag of dependent variable serves as instrument for first lag)
# # Do not reject H0 (i.e. no second order serial correlation) if p > 0.05
# 
# cm <- c("lag(logtrade, 1 × 2)1" = "Lag ln(dyadic trade), 1",
#         "lag(logtrade, 1 × 2)2" = "Lag ln(dyadic trade), 2",
#         "lag(weight, 0 × 2)0 " = "ln(TC connections)",
#         "lag(weight, 0 × 2)1" = "Lag ln(TC connections), 1",
#         "lag(weight, 0 × 2)2" = "Lag ln(TC connections), 2",
#         "lag(gdpcap_d, 0 × 1)0" = "GDP per capita (exporter)",
#         "lag(gdpcap_d, 0 × 1)1" = "Lag GDP per capita (exporter), 1",
#         "lag(gdpcap_o, 0 × 1)0" = "GDP per capita (importer)",
#         "lag(gdpcap_o, 0 × 1)1" = "Lag GDP per capita (importer), 1",
#         "lag(log(pop_d), 0 × 1)0" = "ln(population) (exporter)",
#         "lag(log(pop_d), 0 × 1)1" = "Lag ln(population) (exporter), 1",
#         "lag(log(pop_o), 0 × 1)0" = "ln(population) (importer)",
#         "lag(log(pop_o), 0 × 1)1" = "Lag ln(population) (importer), 1",
#         "lag(fta_wto, 0 × 1)0" = "Regional trade agreement",
#         "lag(fta_wto, 0 × 1)1" = "Lag Regional trade agreement, 1",
#         "lag(wto_dyad, 0 × 1)0" = "WTO dyad",
#         "lag(wto_dyad, 0 × 1)1" = "Lag WTO dyad, 1",
#         "lag(dem_dyad, 0 × 1)0" = "Democratic dyad",
#         "lag(dem_dyad, 0 × 1)1" = "Lag democratic dyad, 1",
#         "lag(pta, 0 × 1)0" = "Preferential trade agreements",
#         "lag(pta, 0 × 1)1" = "Lag preferential trade agreements, 1",
#         "lag(comcur, 0 × 1)0" = "Common currency",
#         "lag(comcur, 0 × 1)1" = "Lag common currency, 1",
#         "lag(alliance, 0 × 1)0" = "Alliance",
#         "lag(alliance, 0 × 1)1" = "Lag alliance, 1",
#         "lag(rivalry, 0 × 1)0" = "Strategic rivalry",
#         "lag(rivalry, 0 × 1)1" = "Lag strategic rivalry, 1")
# 
# models <- list(gmm_all0, gmm_all1, gmm_all2, gmm_all3, gmm_all4)
# names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity++", "Gravity+++")
# 
# rows <- tibble::tribble(~start, ~gmm_all0, ~gmm_all1, ~gmm_all2, ~gmm_all3, ~gmm_all4, 
#                         "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
#                         "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011")
# 
# tabsummary <- modelsummary(models, 
#                            fmt = 3,
#                            #coef_map = cm,
#                            #output = "latex",
#                            stars = TRUE,
#                            gof_map = "none",
#                            notes = c("Fixed effects: Dyad and year.", 
#                                      "Clustered standard errors by dyad and year.",
#                                      "Model: Generalized Methods of Moments (GMM)"),
#                            add_rows = rows)
# 
# gmm_summary <- tabsummary %>%
#   add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
#   kable_styling(font_size = 10, full_width = FALSE)
